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THE NASA AMES RESEARCH CENTER ONE- AND TWO-DIMENSIONAL STRATOSPHERIC MODELS 


PART II: THE TWO-DIMENSIONAL MODEL 

R. C. Whitten,* W. J. Borucki,* V. R. Watson,* T. Shimazaki,* 
H. T. Woodward,* C. A. Riegeljt L. A. Capone, t and T. Becker^ 


SUMMARY 


The Ames two-dimensional model of stratospheric trace constituents is 
presented in detail. The derivation of pertinent transport parameters and the 
numerical solution of the species continuity equations, including a technique 
for treating the "stiff" differential equations that represent the chemical 
kinetic terms, and appropriate methods for simulating the diurnal variations 
of the solar zenith angle and species concentrations are discussed. Predicted 
distributions of tracer constituents (ozone, carbon 14, nitric acid) are com- 
pared with observed distributions. 


INTRODUCTION 


In this paper we discuss the two-dimensional atmospheric model developed 
at Ames Research Center: its geometry, techniques for solving the , associated 

photochemical continuity equations, and the derivation and inclusion of 
appropriate transport parameters. Two-dimensional models are justified 
because: (1) zonal spreading and mixing of trace constituents is much more 

rapid than is meridional spreading and mixing, thus validating the assumption 
of zonal symmetry for most applications; (2) meridional bulk velocity and eddy 
diffuslvlty fields can be constructed by choosing them such that distributions 
of several tracers with different "chemistries" are satisfactorily simulated 
while the fields themselves are in rough agreement with observed winds and 
their variances; and (3) reasonably complete sets of chemical species and 
reactions can be accommodated without the excessive demands on computer time 
that would be required by three-dimensional models. Because this article is 
intended to emphasize the mechanics and numerical methods used rather than 
predictions of trace gas distributions and the effects of perturbations 
thereof on the Earth's ozone shield, we restrict our discussion accordingly. 

To complement this review we refer the reader to our recent publication 
(ref. 1) on the effects of supe^rsonlc transport emissions on the ozone layer; 
additional papers describing results for trace gas distributions in the 


*Ames Research Center 

^Dept. of Meteorology, San Jose State Univ. 
t Informatics , Inc. 



ambient atmosphere and for effects on atmospheric ozone of chlorofluorometh- 
anes and of HCl emitted by space shuttle launch vehicles are now being pre- 
pared. We emphasize that our two-dimensional model is under continuing 
development and has not yet reached the same state of maturity as the one- 
dimensional model of Turco and Whitten (ref. 2) . The geometry of the model is 
shown in figure 1. 


80° N 

r'-^okm 



80° s 


25 VERTICAL GRID POINTS 
AT 2.5 km INTERVALS 
33 MERIDIONAL GRID POINTS 
AT 5° INTERVALS 


The chemistry employed in this 
model is a somewhat simplified ver- 
sion of that used in our one- 
dimensional model studies (e.g., 
ref. 2) . Because of the computa- 
tional time demands of a two- 
dimensional model, it is necessary 
to include only that set of species 
and chemical reactions that is 
essential to realistic simulation of 
species distributions. Hence, our 
chemical set is smaller than that 
in the one-dimensional model 
(ref. 2). In arriving at the chem- 
ical reaction set in the two- 
dimensional model we used simula- 


tions from the one-dimensional model 
Figure 1.- Model geometry. as a guide, ignoring those reactions 

that have very little effect on pro- 
duction and loss rates of the constituents. At each .time step, the model cal- 
culates the concentrations of Oj, 0(^P), 0(^D), N, N02> NO, NOj, N2O, N^Og, 
HNO3, HNO2, H2O2, HO2, OH, H, Cl, CIO, HCl, CFCI3, and CF2CI2. The profiles 
of N2, O2 , CO2, CO, CHi^ , and H2O are held fixed (at experimentally determined 
values) during the calculation. . However, the program is now being modified to 
allow CHj^ and CO, and stratospheric H2O to vary during the computations. 


SPECIES CONTINUITY EQUATIONS 


The model is based on the governing equations for trace constituents: 


$. = P. - L. 


with the flux, $.£ , expressed as the sum of a term that represents large scale 
mean meridional bulk motion and a term that represents smaller scale mixing 
caused by "eddy" motions 



Here, and are the concentration and flux, respectively, of. the ith. con- 
stituent; and are the respective chemical production and loss terms; 
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H and T are the atmospheric scale height and temperature; is a vertical 
unit vector; v is the atmospheric bulk velocity; and Kg is the eddy diffusiv- 
ity tensor for which the symmetry is assumed. The quantities E and T are 
related by the equality 


E = 



(3) 


in which k-Q is Boltzmann's constant, g is gravitational acceleration, and m 
is the mean molecular mass in the well-mixed atmosphere. The temperature, T, 
which is a function of season and latitude as well as altitude, is taken from 
Oort and Rasmussen (ref. 3) for altitudes below 20 km and from Nastrom and 
Belmont (ref. 4) for higher altitudes. The two-dimensional "del" operator is 
represented in the meridional plane of a sphere by 


V 


R cos 6 36 


cos 


6 + 


3s 


(4) 


where 6 is the latitude, bq is a tangential unit vector, and R is the radius 
of Earth. 


The term - h£, which represents the net production rate of a species, 
can be written 





L. 


= E 


e . .^k .,n .n^ 
^Jk ok 0 k 


(5) 


where is ±1, depending on whether the species is produced or destroyed, 

respectively, and k^j^ is the rate coefficient for the reaction. If the 
process is photolytxc, we have nj = 1 and 



( 6 ) 


with J-j^ the photolysis rate. 

Transport parameters are specified in the model for each season. They 
are of two types: bulk velocity, which represents large scale mean meridional 

motion, and "eddy dif fuslvity ," which represents smaller scale motions. The 
parameters that we have adopted are physically reasonable and lead to satis- 
factory predictions of observed tracer distributions. 

The mean meridional circulation is obtained by the kinematic method from 
the averaged equation of mass continuity, which states that 

1^ + V • (pk) = 0 (7) 

dt w 
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where p is the atmospheric bulk density and is the large scale atmospheric 
velocity. With the assumption that the density field is in a steady state, 
the approximate form of this equation in spherical coordinates is 

S els 6 ^ ° 

where the overbar denotes an average with respect to time and longitude, and 
V and w are the mean meridional and vertical velocity components. Equa- 
tion (8) implies the existence of a "stream function" for the total flux 
such that 

2irBpv cos 0 = - ; 2nRpij cos 6 = -g (9) 

oS n aW 

If the distributions of p and V are known, the first of equations (9) can be 
integrated vertically to obtain i)j, and M can then be obtained from the second. 
The required boundary conditions on i|) are discussed below. 

The density was obtained from the hydrostatic equation using temperature 
data presented by Oort and Rasmussen (ref. 3) for altitudes to 20 km; for 
higher altitudes mean (rocket) temperatures were used (ref. 4). 

The V components were based on the observed data of Oort and Rasmussen 
(ref. 3) for altitudes to 20 km. These values were then extrapolated to the 
top of the model by imposition of a simple vertical profile that matched the 
20-km value and satisfied the kinematic constraint that there be no net mass 
flux across any vertical latitude wall. In addition, a meridional behavior 
was imposed such that a three-cell circulation structure (80“ N to 80° S) in 
the stratosphere and lower mesosphere resulted during the summer and winter 
seasons, and a four-cell structure in spring and fall. This should approxi- 
mate the circulation patterns obtained by other investigators (e.g. , 

Murgatroyd and Singleton, ref. 5; Vincent, ref. 6). 

The first of equations (9) was then integrated vertically with the condi- 
tion that = 0 on all boundaries. Meridional wind speeds are generally less 
than 1 m sec“^, except in the upper equatorial troposphere. Small mean 
meridional velocities are, of course, expected because they result from the 
(small) ageostrophic part of the circulation system. Vertical wind speeds 
attain maxima in excess of 4 mm sec~^ in the equatorial troposphere, but are 
usually less than 1 mm sec“^ except in the vicinity of cell boundaries. 

Stream functions for fall in the Northern Hemisphere are shown in figure 2. 

The eddy fluxes are modeled by diffusion coefficients and 

(the local "y" axis is tangent to the meridian; i.e., dy = Bdd) . Because 
the motivation for choosing the general shape of the profile is described 

in great detail in Turco and Whitten (ref. 2), we do not repeat it here. It 
is sufficient to note that the magnitude of in this model is generally 

somewhat smaller than that in a one-dimensional model because vertical trans- 
port in two-dimensional models is affected by large scale bulk motion (w) as 
well as by eddy mixing. The values of B,.,. , which is the horizontal analogue 
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of Kzzi strongly latitude- 

dependent, but height- independent , 
as shown in figure 3; we note that it 
was not necessary to make Kyy height- 
dependent in order to match observed 
tracer distributions. The magnitude 
tends to be large at high latitudes 
in winter and spring when the wind 
shears arising from the polar vortex 
cause the formation of strong eddy 
motions. On the other hand, Kyy is 
smaller in the tropics where zonal 
flow and thus eddy formation are 
weaker. Our meridional variation of 


Kyy is different from that of Louis 
(ref. 7), who adopted larger Kyy at 
low latitude than at high latitudes; 


Louis's choice of Kyy does not yield 
the observed ozone distribution when 


chemistry is included in the model 
(see the results section) . 


The most interesting of the 
however, is the off-diagonal Ky^ 
because it leads to an apparent 
counter-gradient flux of ozone. It 
has been convincingly argued by 
Newell (ref. 8) that the "mixing sur- 
faces" in the lower stratosphere 
slope more steeply toward the pole 
than do the surfaces of potential tem- 
perature. That is, parcels of ozone 
moving poleward in the lower strato- 
sphere tend also to move farther down- 
ward in the transport region to alti- 
tudes where their destruction becomes 



Figure 2.- Stream function; fall in 
the Northern Hemisphere. 



less efficient. Hence we see a 

buildup of ozone column density at Figure 3.- ^yy\ fall in the Northern 

high latitudes. The reason for the Hemisphere. Kyy is independent of 

difference in slope of the potential height in our model, 

temperature and mixing surfaces lies 

in the effect of eddy motions in the lower atmosphere. In poleward moving air 
these motions tend to force the air parcels downward below the surfaces of 
potential temperature, the surfaces of natural movement in quiet stable air. 
One can in principle calculate the K^j by computing the variances of observed 
winds over a long period of time. In practice, one can obtain only order-of- 
magnitude estimates in this way. Our K{_j, which are based on observed tracer 
distributions, are consistent with the order-of-magnitude estimates obtained 


5 





from wind variance calculations. In order to facilitate computation of the 
K-ij, we employ a suggestion by Brasseur^ that they be expressed in the form 

K. . = Q. .(0)Z. .(z) (10) 

-z-J 

where 0 is a function of latitude alone and Z is a function of height alone 
(note that our Zyy = 1). These functions are readily adjusted to yield 
Z-fields that result in smooth distributions of trace constituents. 


BOUNDARY CONDITIONS 


End boundaries are taken at 80° S and 80° N because meridional fluxes are 
expected to be small at those latitudes. Hence, the end boundary conditions 
are taken to be zero flux of all constituents across the vertical boundaries. 
The upper boundary conditions are given by setting the fluxes equal to zero 
for all species. We have checked the upper boundary conditions (at 60 km) for 
the various species against results from a one-dimensional model that extended 
up to 120 km. It was found that the effect of choosing a mixing equilibrium 
condition at the upper boundary had very little effect on any of the constitu- 
ents at altitudes below 50 km. The lower boundary condition used for all 
species except N2O, CH^, HNO3, NO2, O3, HCl, H2O2 and the halocarbons, is 
chemical equilibrium because of their short lifetimes against chemical loss. 
Because HNO3, HCl, and H2O2 are water soluble, their number densities are set 
equal to zero at the lower boundary. The number densities of CHt^, NO2, and 
N2O are fixed at 3.7x10^^ cm”^, 3x10^ cm“^, and 7.5x10^^ cm“^, respectively, 
at the lower boundary, and that of O3 is fixed at 6x10^^ cm“^ (ref. 9) in 
order to conform to measured values. The boundary conditions for CFCI3 and 
CF2CI2 time-variable fluxes that are set equal to the chlorofluorocarbon 
release rates (ref. 10). 


NUMERICAL SOLUTION OF THE SPECIES CONTINUITY EQUATIONS 


We now turn our attention to the numerical solution of the governing 
equations. The method of "time splitting" is used to facilitate the computa- 
tion; that is, the chemistry and transport portions of the equations are in 
general solved at different time steps. The technique for the "chemistry" 
step is discussed first; the technique that applies to transport is described 
in following sections. 


^Private communication from Brasseur, Institut d'Aeronomie Spatiale de 
Belgique, Brussels. 
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We avoid the instabilities associated with "stiffness"^ by employing an 
approach that is closely related to the "families" method developed by Turco 
and Whitten (refs. 2, 11). We prefer the following technique to that of Turco 
and Whitten because in a two-dimensional model it appears to cause fewer dif- 
ficulties at high latitudes where the solar illumination vanishes near the 
winter pole. 


The chemical rate equations are solved using the following type of 
implicit technique. Equations (1) and (5) yield the finite difference form 



At 



( 11 ) 


in which the superscript j indicates the index of the time step At. Equa- 
tion (11) can .be linearized by taking the first term in a Taylor series expan- 
sion about : 


rP.^^ - nP 

-i_ 

At 





( 12 ) 


The m members of the set of mass conservation equations are thus coupled and 

require solution by inversion of a large matrix at each grid point. The 

24 species in the Ames two-dimensional model are grouped into sets of closely 
coupled species that interchange rapidly (small time constants); but the 
species in a set are only weakly coupled to the species not contained in the 

set (large time constants) . For example, a set containing the odd-oxygen 

species O3, 0(^D), and 0 need not contain NO3, N2OJ3, HNO3, HNO2, N2O, or H2O2 . 
Because the computation time required for matrix inversion is roughly propor- 
tional to the square of the number of species, a careful segmentation into a 
number of smaller matrices substantially reduces the computer time required 
for solution. This method allows large time steps (about 1 day for most 
applications) to be taken even though the equations are very stiff. The tech- 
nique is very efficient and yields results that agree well with those obtained 
with a linearized fully implicit method (ref. 1). Table 1 lists the various 
closely-coupled sets and their members. 


The photodissociation coefficient of the ith constituent at wavelength 
A and height z is calculated from the expression 

J.(a,A) = £.(X)o.(A)J„(A)e‘''^^’^^ (13) 


system of "stiff" coupled differential equations is characterized by 
time constants that span a large range of values. Use of large time steps 
does not permit the systems to respond to the small time constants. This 
characteristic can lead to computational instability unless special precau- 
tions, such as that described here, are taken. 
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TABLE 1.- THE SETS OF SPECIES 


No. 


Components 

1 

O 3 , 0(3p) 

, 0(1d), NO 2 , Cl, CIO, HCl 

2 

OH, HO 2 , 

H 2 O 2 , H, HNO 3 , HNO 2 

3 

NO, NO 2 , 

NO 3 , N, N 2 O 5 , HNO 3 , HNO 2 


Note: Solution of the matrix equations must 

be in the order listed; for example, 
NOx compounds (set no. 3) must be 
determined after the solution for set 
no. 2 , not no. 1 . 


where is the absorption cross section of the ith constituent, is the 
intensity of solar radiation (number of photons cm~^ sec"^ (unit A)”^ at the 
top of the atmosphere), e^(A) is the quantum yield, and x is the optical depth 
given approximately by 


t(A,3) 






ds + a (A) 



ds + o 


NO, 


(A) 



“I 


ds 


sec 


X 

(14) 


for small solar zenith angle, x* At larger x a tay tracing routine is used to 
determine the physical path length. The optical depth is, of course, implic- 
itly time-dependent through temporal changes in ^ 03 ’ clarity we have 

omitted time (t) from the arguments. In order to calculate the integrated 
dissociation coefficient over a wavelength region (Aq to Ay) , we divide it 
into a number of intervals (A^ “ ~ ^ 2 ’ * • > which include 

the Schumann- Run ge bands. Then the total photodissociation coefficient is 
calculated by 


J.(B) 



J^(A,s)dA 



-^oo(^)o • (^)exp 

O' O' 



(A,s) + T 

*^3 


(A,s) + T 


NO, 



dA 

(15) 


Because Oq^ and hence Tq are very strong functions of A in the Schumann- 
Runge band system (175-205 nm) , special methods must be used in that spectral 
region. Hudson and Mahle (ref. 12) performed a line-by-line integration using 
theoretical line profiles for column densities N between 1x10^^ cm ^ and 
7x10^^ cm“^ and for temperatures between 150 and 300 K. Interpolation con- 
stants were then found for 19 wavelength bands spanning the 175-205 nm spec- 
tral region for both the oxygen photolysis rate and the optical depth. 

Because the 5538 interpolation constants must be stored in memory, and because 
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interpolation is required to obtain a specific value of the photodissociation 
rate and optical depth, the Hudson-Mahle model is expensive to use. 

In order to reduce the computational time Shimazaki and Ogawa (refs. 13 
and 14) developed a modification of the Hudson and Mahle technique in which 
they derived polynomial formulas of the seventh degree for computing the inte- 
grals, which appear in equation (14), as functions of N: 

1 

£n I = P, (£n iV)^ (16) 

fe =0 ^ 

The coefficients P^ are determined by a least squares best fit to the Hudson- 
Mahle results. The Ames two-dimensional model uses 6 equal intervals of 5 nm 
each in the Schumann-Runge bands and 51 intervals throughout the remainder of 
the solar spectrum. Shimazaki et al. (ref. 15) have carefully compared the 
obtained in this way for all of the important stratospheric photodissocia- 
tion processes, and found excellent agreement (within 10 percent) with the 
Hudson-Mahle approach. 

Ideally, one would like to perform computations of the concentrations of 
stratospheric trace constituents that fully include the effects of diurnal 
variations. Unfortunately, the small time-steps required for a full diurnal 
simulation are as a rule much too expensive in terms of computer time require- 
ments. Hence, we prefer to use some means of approximating the diurnal 
effects. Such a simulation will allow for change of the solar zenith angle 
during the day, absence of solar illumination during the night, and very sub- 
stantial change in concentration of some of the constituents from day to 
night. The diurnally-averaged photodissociation rates are obtained by approx- 
imating the time integral of equation (15) (see Cogley and Borucki, ref. 16). 
The approximation considers the number of hours of daylight at each latitude 
and season, and matches various properties of the exponential integral and 
exponential functions . 

Turco and Whitten (refs. 2, 17) have developed a method for averaging the 
species concentrations which simulates quite well the results obtained from a 
full diurnal calculation of the trace constituents. The most sensitive test 
of the technique is the calculation of the N 2 O 5 abundance because of the slow 
build-up of N 2 O 5 over many diurnal cycles. The results obtained with the 
Turco and Whitten method are within 20 percent of those obtained with a full 
diurnal calculation. In applying this averaging technique to the two- 
dimensional model, one must be sure to use it at evevy latitude because of the 
latitudinal dependence, particularly at high latitude, of day-night ratios. 

Solutions for the transport processes are obtained from a mass-conserving, 
forward-time, space-centered finite difference formulation of the equations. 

The transport computations are further time split into those for vertical and 
horizontal advection and diffusion. A backward-time implicit formulation is 
used for the horizontal diffusion and a forward explicit formulation is used 
for all other transport processes. The calculations for advection are modi- 
fied by a term M that improves the accuracy from first order to nearly second 
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order. This technique for improving the accuracy is practical because the 
velocities are prescribed and do not vary daily. 


The differential and corresponding finite difference forms of the hori- 
zontal transport equations are as follows. 

Advect-ion (forward-time, space-centered, explicit) : 

(17a) 
(17b) 

where 


U 




. - 


cos 0 




,Z^ 

> Q "2'“^“ 1 


- V 


■Z^ 

t-\,g ^-l 


Z^ . = Z(6., 2 ., t, ) 

V . . = cos 6 . y (6 . , 2 .) 


Diffusion (forward-time, space-centered, implicit): 


JZ ^ 1 3 

3t cos 6 i?36 


^os oj. 


Z* + K 


yy i?80 i/2L92 


rp+(4 


e-iSHi) 


At 


y^+1 s 7^ + 

^ig' ig cos 0. 2(i?A0)2 

't' 


'K* 

yy 

L\i+i ,g 


+ z* Vz^'!'^ . - 


yy 

ig/ 


^+1 




- fz* + z* 
yy yy 
ig d-\,g! 


- z^+' . 


1 / At \ ~ // 

cos 0 . \4ZA0A2 j yz 1 i+\ ,g- 


- Z^ 


1 At 
cos 0. 2ZA0 


- Z^ ^ 

,j+l ^+l,J-li 


z 


y^ 


y t-1 , 


j+i 


(18a) 


,j^+i ,J %z^Th-l ,J 


(18b) 


where 


z* 

= cos 

0 . 

Z* 1 

yy 



yy 

K* 

= K 

+ 

M 

yy 

yy 



k 

= cos 

0 . 

K 

yz 



yz 
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T . - T 

1 . 1 

H. . T. . Ikz 


M = 


At 


max 


u? . ,y? .,y? 


. ,y? . ,y? 

j+i ^,j-l 


y? . ,y? . ,y 


2 .,y2 . I 

y+i,j’ ^+l,J+lJ 


Equations (18b) can be rewritten as a tridiagonal system and solved for 


by a standard inversion technique (e.g., Ames, ref. 18; also see Turco and 
Whitten, ref. 2). 


The differential and corresponding finite difference forms of the verti- 
cal transport equations are as follows. 

Advection (forward-time , space-centered, explicit): 


U 3 , 


"Z-.J V^Vl "Z-.J+l "Z-jJ-l ^-..7-lJ 

Diffusion (forward-time, space-centered, explicit): 

3^ 


3t 3a 1 zz 


L FM + /l + 1 1 1 ^ 

I a2[32 yff T da/'^J ^ ya 36 
At 


^,J 2(Aa)^ 


I aa 

' 7^ . .1 ■ 


+ t: 


L> ^,J+l 


aa 

i,J/ 


W + I* . 


,i) 


I aa 


+ X Va^ . + . 

aa II 

i,0-l'^ 


,J-1^ 


At 

2Aa 


^^zz^PiJ+l^iyj+l 


^^zz^T^i 


,J-1 


At 


Ai?A6Aa 


ya lt+l,j+i t-i,j+i 


- K . 

ya I ^+l ,j 


- 

-1 


(19a) 

(19b) 


(20a) 


(20b) 


In our model we have chosen A6 = 5° and Aa = 2.5 km. Because the formulas for 
the finite difference forms of the boundary condition are very lengthy, but 
not essential for immediate considerations, we present them in the appendix. 
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It is important to note that our differencing techniques for transport 
are mass-conserving. They are also quite accurate as we have verified with 
the aid of simple sinusoidal functions. ■ There proved to be very little dis- 
tortion or change in amplitude of the sinusoidal forms as a result of advec- 
tive effects. A form of X with sinusoidal time and space dependence was 
inserted into the finite difference forms of the advective transport equa- 
tions. Because the differential forms predict that such waves will not be 
distorted nor their amplitudes changed with time, we desire the same type of 
behavior for the finite difference solutions. Wave distortion and attenuation 
proved to be small. 

Turco and Whitten (ref. 17) showed that in the stratosphere, diurnal 
averaging of the vertical transport terms in the continuity equations can be 
neglected. This is so because the time copstants of species that experience 
large day-night variations are small compared to typical transport times. It 
is easily shown that similar considerations hold for the horizontal and "off- 
diagonal" transport. Hence we can neglect diurnal averaging in computing 
transport of species. 


SOME RESULTS 


The Ames two-dimensional model was used to study the distribution of 
various trace species (ozone, nitric acid, and carbon 14) in the stratosphere 
and to assess the effect on ozone of various pollutants. In this section we 
compare our predictions of ozone, nitric acid, and carbon 14 with measurement 
and show how our set of transport parameters is supported by observation. 

Figure 4 shows some calculated profiles of ozone concentration for the 


midspring 

season together 

with measured 


WILCOX et al. (ref. I9) 
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Figure 4.- Vertical distributions of 
ozone at various latitudes. 


data summarized by Wilcox et al. 

(ref. 19) for altitudes below 35 km; 
the high altitude data were obtained 
by polar-orbiting satellites 
(I. Eberstein, Goddard Space Flight 
Center, private communication). The 
topside concentrations (above 
40 km) , which are quite close to the 
measured values, are strongly depen- 
dent on the rate of ozone destruc- 
tion by HOjj catalysis; OH concentra- 
tions in excess of 10^ cm~^ are 
required to match the observed ozone 
profiles. At altitudes below 30 km, 
where ozone abundance is dominated 
by transport, our predicted profiles 
are quite close to those observed. 

In comparing the predicted 
ozone abundance in the transport 
region with observed data it is 
instructive to employ the vertical 
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column density as a function of latitude. Transport parameters must, of 
course, be selected such that the predicted column densities are reasonably 
close to observed values. Louis (ref. 7) attempted to do this with the aid of 
a time-dependent model whose eddy diffusion coefficients were deduced from 
observed ozone distributions. Then his model was tested by simulating the 
distribution of radioactive debris from Chinese and American nuclear explo- 
sions. As one might expect, the agreement was very good because chemical pro- 
duction and loss of ozone was not included in the calculations used to deduce 
the K's from observed ozone distributions. In other words, ozone, like the 
radioactive debris, was treated as an inert tracer. 


In order to check the effects 
when ozone chemistry is included, we 
have run our model (which includes 
photochemical production and loss 
terms) using Louis's transport param- 
eters. The resulting ozone column 
density distribution (as a function of 
latitude), shown in figure 5, is 
obviously in strong disagreement with 
observation. However, when the same 
type of distribution is calculated 
with our set of transport parameters, 
the resulting column density distri- 
bution, shown in figure 6, is in very 
good agreement with observation. This 
emphasizes the need to take into 
account the complete chemistry, 
sources, and sinks of any tracer used 
to deduce transport parameters. 



(a) Summer in the Northern Hemisphere. 


CM 



Figure 5.- Ozone column density as a 
function of latitude, using Louis' 
transport parameters. 



(b) Fall in the Northern Hemisphere. 


Figure 6.- Ozone column density as a function of latitude; using our 

transport parameters. 
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Our predicted nitric acid column distribution for the same seasons is 
illustrated in figure 7. Since the chemistry of NO^ compounds differs from 
ozone chemistry, we expect the distributions to be different. In particular 
we expect HNO 3 column density latitudinal gradients to be larger than the 
corresponding ozone gradients because HNO^ loss rates at all latitudes are 
smaller than those of ozone. This difference is both predicted and observed. 


Finally, figure 8 compares the predicted and measured carbon 14 concen- 
trations. Although the removal rate of from the lower stratosphere 
appears to be somewhat too rapid, insistence on better agreement is probably 
not reasonable because the quality of the observed distributions at altitudes 
above 20 km is poor. We believe that we have found quite good overall agree- 
ment with the set of tracers (ozone, HNO 3 , and carbon 14). 


MURCRAYetal. (ref. 2i) 

NO. HEMISPHERE BALLOON 


SPRING 

□ 

JAN 

12-16 

■ 

MAY 

12 

WINTER 

0 
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18 

• 
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12 

SUMMER 

A 

APR 

15 

♦ 

SEP 

12 


0 

APR 

18 

A 

NOV 

12 



Figure 7.- Nitric acid column density 
as a function of latitude. Balloon 
measurements were made in 1970. 


I YEAR RUN 


PREDICTIONS 


MEASUREMENTS 


0 °N. LAT 

30 

60 


□ 0 °N. LAT 

o 30 

A 60 



Figure 8 .- Vertical distributions of 
carbon 14 at various latitudes. 
Experimental points are from 
Johnston et al. (ref. 22). 


CONCLUSIONS 


We have presented a complete and concise outline of our two-dimensional 
atmospheric computer simulation. Some of the unique aspects of our model are: 

1. We are able to match simultaneously the observed distributions of 
several stratospheric tracers. 

2. Our transport parameters are physically realistic. 

3. We use an averaging scheme for simulating the effects of diurnal 
variations of atmospheric constituent concentrations. 
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4. We use an accurate computational technique for transport of species. 

5. Our model is computationally efficient. 


Ames Research Center 

National Aeronautics and Space Administration 

Moffett Field, California 94035, April 14, 1977 
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APPENDIX 


BOUNDARY CONDITIONS 


The finite difference forms of the transport equations for the boundary 
points are as follows. 


Horizontal Advection 


Boundary at 80° S — zero flux boundary condition: 


^+1 = + At 

1 , j l,j 2b.y cos 6 


- V .+ y .1^ ~\ 

1 (_ 2,J 2,J 1 1 


Boundary at 80° N — zero flux boundary condition: 


^+1 ht 

N,j Nj 2Az/ cos 0 


[V 1 ,/ n - 1 , j , j] 


Horizontal Diffusion 

Boundary at 80° S — zero flux boundary condition: 

1,J 1,J 2(Az/)2 cos 0 I yy yy 2,j l ,j 

V2,J 1,J7\ 


hb,yhz cos 0 


L2,j' ' l,j' 'J 




Boundary at 80° N — zero flux boundary condition; 


. (y ^ y 

2(Ai/) 2 cos 0 [ yy yy 

\N,j N- 


hbiyhz cos 0. 


^ (k + k - 

cos 0 I yy yy Jl N 

\N,j N-\,o/\ 

■ \z fe-l,j+l “ 

~ k 

2Az/ cos 0^, 


• 
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Vertical Advection 


Upper boundary (j = M) — zero flux boundary condition: 


^+1 ht 

i,M iM 2Aa 


+ l/. ^ 

I zM ^M ^M- 1 


Lower boundary (j = 1) — specified density: 


yk+1 




density specified on chemical equilibrium 


value calculated from chemical equations or specified flux 


(The specified flux could be added in the diffusion transport calculation — 
that is, modelled as a diffusion transport, or added in the advection trans- 
port, or added in the advection transport calculation. In this program, the 
specified flux is added in the diffusion transport calculation and the flux 
due to advection transport is set equal to zero.) 


t , 1 t , 1 


As 


L /. . 

I t , 1 it. 



Vertical Diffusion 

Upper boundary (j = M) — zero flux boundary condition: 


/+1 ^ yk _ At 
i,M ^iM 


K \ + K // - I 

yzyi+l,M yz 


\(k + k 

2 (As) 2 y ss ss ly iM %M-l j\ 
iM '-I 

(a: hj + (K hj ... ^ 

j ss T ZM ^M ss T zM-\ zM-\\ 


At 

2As 


Lower boundary (j = 1) — specified density: 


yk+l 

^,1 


yk+\ 


density specified or chemical equilibrium. 

value calculated from chemical equations or specified flux. 


17 



hz Ihljhz 


4-1,2. 


^l 


^,2' 


At 


(As)2 


+ a: |fc , - ^ , 

l 22 22 /I ^ ,2 ^ , 1 . 

L\ii i2> 


+ 


A2 


(is: iz„) . ^ + (X ?z„) . 

zz T zl z,i zz T' zi ^,; 


where s|)g = flux specified at the lower boundary. 
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